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Systems of fermions described by the three-dimensional (3D) repulsive Hubbard model on a cubic 
lattice have recently attracted considerable attention due to their possible experimental realiza- 
tion via cold atoms in an optical lattice. Because analytical and numerical results are limited 
away from half-filling, we study the ground state of the doped system from weak to intermediate 
interaction strengths within the generalized Hartree-Fock approximation. The exact solution to 
the self-consistent-field equations in the thermodynamic limit is obtained and the ground state is 
shown to exhibit antiferromagnetic order and incommensurate spin-density waves (SDW). At low 
interaction strengths, the SDW has unidirectional character with a leading wave-vector along the 
(lOO)-direction, and the system is metallic. As the interaction increases, the system undergoes a 
simultaneous structural and metal-to-insulator transition to a unidirectional SDW state along the 
(11 Indirection, with a different wavelength. We systematically determine the real- and momentum- 
space properties of these states. The crossover from 3D to two-dimensions (2D) is then studied 
by varying the inter-plane hopping amplitude, which can be experimentally realized by tuning the 
distance between a stack of square-lattice layers. Detailed comparisons are made between the exact 
numerical results and predictions from the pairing model, a variational ansatz based on the pair- 
ing of spins in the vicinity of the Fermi surface. Most of the numerical results can be understood 
quantitatively from the ansatz, which provides a simple picture for the nature of the SDW states. 

PACS numbers: 75.30.Fv, 71.15.Ap, 71.45.Lr, 71.10.Fd, 75.10.Lp, 75.50.Ee 



I. INTRODUCTION 

Over the past several years, optical lattices have be- 
come an increasingly powerful tool for emulating many 
systems in condensed matter physics^^ An optical lat- 
tice can provide exceptionally clean access to a vari- 
ety of model many-body Hamiltonians in which param- 
eters can be systematically tuned and controlled. Thus, 
they make possible quantitative experimental study of 
the properties of interacting electron models, which have 
proven extremely challenging for analytic and numerical 
approaches alone. The combination of these approaches 
presents unprecedented opportunities for improving our 
understanding of interacting electron systems, by test- 
ing theoretical concepts and increasing the accuracy and 
predictive power of numerical approaches via comparison 
with experiment. 

The one-band Hubbard model is one of the most fun- 
damental models in condensed matter physics. It has 
been widely studied in two dimensions (2D]P^ as the 
simplest model for the Cu-0 plane in cuprate supercon- 
ductors. For the three-dimensional (3D) Hubbard model, 
however, considerably less is known from both theoretical 
and experimental sides. Optical lattices play, in this re- 
spect, a particularly fundamental role as they allow for a 
clean experimental realization of the 3D model and offer 
the interesting possibility of tuning the hopping parame- 
ter along one direction, t±, thereby allowing a systematic 
study of the evolution of properties as the system crosses 
over from 3D to 2D. 

Thanks to advances in the ability to directly cool atoms 
in optical lattices^, experiments are nearing the realiza- 
tion of phases with magnetic order. It is thus particu- 



larly important and timely to understand such phases in 
the Hubbard model. Somewhat surprisingly, apart from 
half- filling (one particle per site), which displays anti- 
ferromagnetic (AFM) order, the nature of the magnetic 
properties in the 3D Hubbard model has not been char- 
acterized, even at the mean- field level. In this work, we 
study the magnetic properties in the ground states of 
the 3D Hubbard model and in the crossover regime, us- 
ing generalized Hartree-Fock (HF) theory. It is shown 
that the system has a tendency to form unidirectional 
spin-density wave (SDW) states with AFM order and a 
modulating wave along either the (100)- (at low U/t) 
or the (11 Indirection (at higher U/t). We examine the 
evolution of the SDW wavelength in the full mean-field 
solution as J7, density and t± vary and characterize the 
ground state by its properties in real and momentum 
space. 

Despite the simple nature of the mean-field approach, 
the determination of the correct eq uilibri um properties in 
these models is not straight for war cP^J The main chal- 
lenge lies in finding an unbiased strategy to determine 
the leading wave- vector (s) characterizing the spatial de- 
pendence of the order parameter. Calculations are per- 
formed in a real space simulation cell and most choices 
of the cell will return solutions that are biased by finite- 
size effects. This is often further complicated by shell 
effects and sensitivity of the solution to the topology of 
the Fermi surface, which often lead to local minimum 
solutions. 

To overcome the difficulties, it is necessary to move to 
larger and larger cells and gain insights from the evolu- 
tion of the corresponding solutions. This line of attack 
has become increasingly possible because of the dramatic 
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increase in computing power and continuous algorithmic 
progress. In the present work, we combine such an ap- 
proach with more targeted searches to obtain the global 
minimum solution of the self-consistent-field (SCF) equa- 
tions in the thermodynamic limit. Furthermore, we show 
how the numerical results can be understood by a varia- 
tional ansatz based on the pairing of spins in the vicinity 
of the Fermi surface. Detailed comparisons are made 
between the direct numerical solutions and the pairing 
ansatz predictions. The excellent agreement helps to pro- 
vide a simple, predictive picture for the properties of the 
SDW states. 

The mean-field approach is often the starting point 
in the study of strongly interacting systems such as the 
Hubbard model. Although the approximations involved 
can lead to significant errors, mean-field theory often 
provides insights into qualitative aspects of the behav- 
ior of many-body systems. Moreover, comparisons with 
quantum Monte Carlo results 19 have showrP^ that, in 
2D, the mean-field solution captures the basic physics of 
SDW states at intermediate interaction strengths, and 
provides a good qualitative (or even quantitative in some 
aspects) description of the magnetic correlations in the 
true ground state. Because it is reasonable to expect a 
similar level of accuracy for the models presently consid- 
ered, we expect that our findings will provide guidance to 
many-body approaches and experimental studies alike. 

We have limited our study to U < 6t, where the mean- 
field approach can be expected to offer useful insight. Be- 
low we will discuss the mean-field predictions and their 
implications (and the caveats) on the tru e ma ny-body 
states drawing from the comparison in 2E) 19 * 20 [ Clearly, 
the form of generalized mean-field theory considered in 
this work will not capture exotic instabilities, such as un- 
conventional pairing order. Indeed, as U increases, it will 
become increasingly inadequate for magnetic properties 
as well. 

The remainder of the paper is organized as follows. In 
Sec. [IlJ we introduce the Hamiltonian, and briefly outline 
some of its basic properties to facilitate the ensuing dis- 
cussion. In Sec. |III[ we summarize the strategies used to 
solve the mean-field equations. Results for the 3D model 



the Hubbard Hamiltonian reads 



are presented in Sec. IV numerical results for the (100)- 
and the (lll)-SDW are followed by a discussion where 
the pairing ansatz is first introduced and then used to 
help understand the numerical findings. The dimensional 
crossover results are then presented in Sec. [V| followed 
again by a discussion based on insights from the pairing 
ansatz. We conclude in Sec. IVII 



II. BACKGROUND 

Given our goal to study both the 3D case and the 
crossover from 3D to 2D, it is most convenient to define 
the 3D Hubbard Hamiltonian as a stack of square-lattice 
planes. We will use r = (x, y) to denote in-plane coor- 
dinates and z to label the planes. With this convention 
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where the operator c\ zo . (c rza ) creates (annihilates) a par- 
ticle with spin a (a =t, I) at site (r, z) and n rza is the cor- 
responding number operator. The hopping amplitude t is 
between nearest neighbor sites within a plane (denoted 
by (it') in the summation), t± is the inter-plane hop- 
ping amplitude between nearest neighbor sites belonging 
to different planes (denoted by (zz f ) in the summation), 
and U > is the on-site interacting strength. Through- 
out this work, energy is quoted in units of t and we set 
t = 1. The Hamiltonian in Eq. ([!]) describes the 3D cubic 
Hubbard model when t± = 1, the crossover between the 
square and cubic lattices when < t± < 1 and a stack of 
decoupled 2D Hubbard planes when t±=0. Only unpo- 
larized systems are considered, i.e., the average densities 
of the two spin species are kept equal: = n±. The 
nature of the ground state is thus characterized by three 
parameters: the inter-plane hopping amplitude t±, the 
on-site repulsion U and the doping (hole density) 



h = l- (n t + 71+). 



(2) 



The particle-hole transformation, c\ za . — ^ (— l) x+yJrZ c rzcn 
maps the h < sector into the h > one, regardless of 
the value of t± or J7, and we therefore confine our study 
to h > 0. 

At half-filling, h = 0, the non-interacting Fermi surface 
is given by —2 (cos k x + cos k y + t± cos k z ) = 0. Despite 
the lack of symmetry between the z- and the r-directions 
for any t± ^ 1, perfect nesting via Q = (±7r, ±7r, ±7r) 
remains throughout the whole surface, and causes an 
AFM instability for any t± ^ and arbitrary small U 
values. The evolution of the non-interacting half-filled 
Fermi surface as t± varies is shown in Fig. [I] In the first 
column, representing the 2D limit, the Fermi surface has 
no dependence on k z and any wave-vector of the form 
(+7T, +7T, q) is perfectly nested on it. The arbitrariness of 
q is reflected in the complete lack of correlation between 
the r-planes. The large nesting degeneracy is abruptly 
interrupted as soon as t± ^ 0, and Q remains the only 
nesting vector as the system evolves from the 2D limit 
toward 3D. The middle and bottom rows illustrate pro- 
jections of the Fermi surface along the (100)- and (111)- 
directions; as we shall see in See's [IV| and \V\ the pro- 
jected area of the Fermi surface plays a central role in 
determining the character of the SDW in the proximity 
offc = 0. 
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are N x N matrices with elements 




FIG. 1: (Color online) Non-interacting half-filled Fermi sur- 
face from different view angles: 3D (top), along [010] (mid- 
dle), along [111] (bottom). From left to right, the columns are 
for t± = 0, 0.5 and 1, respectively. Note that only the surface 
in one octant is shown in the bottom row. Perfect nesting via 
Q = (±7r, ±7r, ±7r) holds for any t± at half- filling. 



III. METHOD 

The following mean-field formalism in real space is 
used in this work. A simulation cell of N sites is defined 
by three vectors, Li, L2 and L3, whose components are 
integers. Bloch states are then introduced as 



c^(k) oc ^2 c /3+l exp [ik • L] , 



(3) 



where L are vectors of the form L = niLi +n<2\j<i + 77,3113, 
k is a reciprocal lattice vector that is free to vary within 
the first Brillouin zone (BZ) defined by the L^'s, and j3 
labels sites inside the simulation cell. Using these states, 
the mean-field Hamiltonian can be decoupled into a sum 
of k-dependent pieces, Hq = 5^ k i?o(k), with each piece 
of the form 



ffo(k) = [c|c 
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where (c^) represents a row of operators c#^(k) 
(c^(k — G)) with index j3 running through the N sites of 
the cell. A non-zero value of G causes the spin densities 
at f3 and /3 + L$ to be related via a rotation by G • L$ 
around the z-axis. Charge and spin densities along z- 
direction obey periodic boundary conditions. H and S 1 * 1 
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(5) 



where £# 7 (k) = exp(zk • L)£/3 j7+ l, and Dp a , Sg and 
\i are determined by the requirement that the free energy 
F = (H)o — TSo is a minimum for the targeted average 
density n = n^ + n±. This amounts to the following SCF 
(gap) equations 
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To locate the ground state we proceed with two com- 
plementary approaches. In the first approach we select 
the L^'s so that they span a large supercell containing 
(9(5000) sites. A twisted boundary conditio n 1 23 * 24 ! is ap- 
plied, namely, using a single randomly selected /c-point 
in place of the integrals in Eq. (J6|. The iterative process 
is started with various initial states, including random 
ones, and multiple annealing cycles are performed. In 
each cycle a random perturbation (whose strength can 
be controlled) is applied to a converged solution and the 
self-consistent process is repeated. Separate calculations 
for different /c-points are done to check for consistency. 

Once an understanding of the character of the ground 
state is gained, we use a second approach to target the 
specific family of states compatible with the results of 
the random search. For instance, suppose the random 
search finds a unidirectional SDW at small U values with 
wave-vector along the (lOO)-direction. We then choose a 
cluster of Li = (L, 0, 0), L 2 = (0, 1, 0) and L 3 = (0, 0, 1) 
with G = 7r((— 1) L+2 ^, 1,1), where I is the number of os- 
cillations of the order parameter, chosen to be an integer 
or half an integer. For a given set of the three parameters 
U , /i), L is finely scanned (with L on the order of 50 
and step size of 1) until the energy minimum is found. 
A large number of /c-points is used (on the order of 100 
in the two short directions and a few in the other) so 
that the character and properties of the targeted states 
can be accurately determined. This approach allows us 
to study different forms of SDW and long wavelength 
modes without increasing the computational cost. 

In our study, we mix the two numerical approaches 
as needed and use them in complementary ways. For 
example, comparison of energies among several families 
of SDW is made with the second approach. To confirm 
the correctness of the ground state, the solutions are then 
checked against different initial states and annealing pro- 
cedures using the first approach on a supercell commen- 
surate with the optimal wavelength. 

Various observables are computed to characterize the 
converged solutions. The local charge density p and the 
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FIG. 2: (Color online) Charge density p (left) and order pa- 
rameter m (right) of the solution for the 3D Hubbard model. 
Shown is a 16 x 16 x 16 supercell, with h = 13/128 at U = 2.5. 
A linear wave is seen along the ^-direction, with uniform AFM 
order in the xy-p\ane. The bottom panel shows a line cut 
along z-direction, with dots the actual data and the line a 
sinusoidal fit. 
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FIG. 3: (Color online) Spatial dependence of the order 
parameter m for h — 13/128, same as in Fig. [2] but with 
U = 2.9, on a 16 x 16 x 16 supercell (left) and a 16 x 16 x 14 
supercell (right). Uniform AFM order in the xy-plane disap- 
pears for L z = 16, but linear SDW along z-direction is seen 
again on the right with L z = 14. 

with oc J] R exp(— ik • R)cR a . We use to compare 
the converged mean-field solution with the pairing ansatz 
prediction and to characterize the Fermi surface of the 
ordered phase. 

IV. 3D RESULTS 
A. SDW correlation in the (lOO)-direction 



local order parameter, identified as the local staggered 
magnetization m, are defined as 



p(R) 
m(R) 



(nrzl)) , 



(7) 
(8) 



and used to characterize the state in real space (here 
R = (r, z)). Since all the minimum energy solutions 
we find are unidirectional spin/charge density waves 
(SDW/CDW), it is natural to characterize them by their 
modulation wavelength along a relevant Cartesian axis, 
Asdw/cdw? denned by the leading component of the 
Fourier transform of p and m, respectively. The min- 
ima in the CDW are found to coincide with nodes in 
the SDW, thus 2Acdw = Asdw- Below we will some- 
times discuss our results in terms of a single wavelength 
A = Acdw, which can also be identified as the distance 
between two consecutive nodes of the order parameter. 
When we refer to the direction of the modulation, we will 
use (100) to denote symmetry-equivalent [100]-directions, 
and similarly for [110] and [111]. 

To characterize the system in momentum space, we 
use the momentum distribution and the momentum- 
resolved single-particle spectral function A^(oj) : denned 
as 



A^ a {uj) = ^Im(c ka (u; - H + E 



irj) 
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At half-filling, the existence of perfect nesting allows 
an AFM solution for any U > 0. Away from half- filling, 
perfect nesting ceases to exist and a finite critical value 
of the interaction is needed to cause the onset of order. 
The critical value U c depends on ft. Using the first ap- 
proach described in Sec. |III[ we have determined that, 
just above U c , the ground state of the system is an SDW 
with modulation along the (lOO)-direction. Figure [^illus- 
trates the spatial dependence of p and m in a 16 x 16 x 16 
supercell at ft = 13/128 ~ 0.10 and U = 2.5. The SDW 
is characterized by a single wave- vector and A^ 10 o) = 8. 
The amplitude of the SDW is 0.1, roughly thirty times 
larger than that of the CDW. The simple form of the 
order found for m(R) is indicative of the proximity of U 
to U c . All of the observations above are consistent with 



the pairing model, as discussed below in Sec. |IV C 
We next examine the evolution of A 



(100) 



as the inter- 



action strength changes. Keeping ft and the simulation 
cell unchanged and increasing U from 2.5 to 2.9, the ran- 
dom search returns the state displayed on the left panel 
of Fig. [3j which suggests that a more complicated type 
of order is seemingly settling in. We apply our second 
approach, using a dense /c-point grid and a 1 x 1 x L sim- 
ulation cell, to search for the optimal wavelength. We use 
large L (containing about 8 nodes in the cell) and vary 
its value until an energy minimum is found. Figure [4] 
shows the result of such minimization in terms of A 



the minimum occurs when A 
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7, indicating that the 
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FIG. 4: (Color online) Energy of (100)-SDW (blue) vs. A (100 ) 
for a system of h = 13/128 and U — 2.9. Horizontal lines are 
the energies of the calculations shown in Fig. [3] The minimum 
of (100)-SDW is reached when A(ioo) — 7. The state in the left 
panel of Fig. [3] leads to an energy higher than the minimum 
but lower than the energy with A(ioo) =8. 



16 x 16 x 16 supercell is not commensurate with the wave- 
length of the minimum energy SDW state and that the 
pattern in the left panel of Fig. [3] is a result of frustration 
from an incommensurate supercell size. We next return 
to our first approach, and perform a new mean-field cal- 
culation, with random initial guess and annealing, on a 
16 x 16 x 14 supercell, a size which is commensurate with 
the wavelength of the minimum energy solution. And in- 
deed we find the predicted state with A^ o) = 7 correctly 
reproduced (right panel in Fig. [3|. 

In Fig. |4j we report the energies of the two large su- 
percell calculations of Fig. [3J to verify that the energy 
obtained in the 16 x 16 x 14 supercell is correctly repro- 
duced by the lxlxl cluster search with A( 10 o) = 7. 
The energy of the 16 x 16 x 16 calculation, on the other 
hand, falls between those of A^ 10 o) = 7 and 8. This gives 
a clear illustration of the characteristics of the two types 
of approaches. The supercell being incommensurate pre- 
vents the solution from collapsing onto the lowest energy 



SDW state of A 



(100) 



7. The self-consistent solution in 



a large supercell then finds a different pattern that cor- 
responds to the true ground state compatible with the 
imposed constraint. The energy of this state, computed 
by fixing the density and converging the value using a 
dense /c-point mesh, is higher than the global minimum, 
but lower than that of the SDW state with an imposed 
wavelength A^ 10 o) = 8- 

We proceed to determine the exact dependence of the 
wavelengths on h and U by explicit solutions of the SCF 
equations in 1 x 1 x L clusters. Verifications of the results 
are done on large supercells whose sizes are commensu- 
rate with the wavelengths. Our results are summarized 



FIG. 5: (Color online) Characteristic wavelength as a func- 
tion of U at various doping, a^ioo) gives the modulation wave- 
length, A(ioo), in units of 1/h. As U is increased, the value 
of ck(ioo) converges to approximately 2/3 at small h (slightly 
larger at larger h). 



in Fig. [5j 



with 



^(100) 



defined as 



^(100) = ^A(ioo)- 



(ii) 



When the doping is small, the wavelength of the modula- 
tion is proportional to 1/h, with a^ioo) almost indepen- 
dent of U and roughly equal to 2/3. For larger h, a(ioo) 
converges to a slightly larger value. There is a general 
trend of an increase of c^ 10 o) as U approaches U c from 
above. We will be able to rationalize these trends within 



the pairing model in Sec. IV C 



The evolution of the properties of the 3D SDW state 
with interaction U is similar to what is observed in 2D. 
Figure [6] shows ID cuts of m and p in the z-direction, at 



h = 0.05 and U values such that a. 



100) 



has saturated to 



~ 2/3. The figure illustrates the existence of U c and the 
increase of the SDW and CDW amplitudes with U. It 
also shows the crossover from a regime where the order 
parameter has a smooth sinusoidal modulation and the 
holes are delocalized, to one where it is characterized by 
domain walls or stripes, with holes localized in the nodal 
regions. 

There exists an important difference in the physics of 
the mean-field ground state of the 3D system and its 2D 
counterpart. While, in the latter, the system remains in- 
sulating when lightly doped, the 3D model immediately 
turns metallic. The difference is a consequence of the 
different behaviors of the modulating wavelength. In 2D, 
a is unity, independent of h and U, while in 3D it varies 
with parameters and has a non-integer value. To illus- 
trate this in a simple case, consider a value of doping h 
such that A (= a/h) is an integer. For such a system to 
be an insulator, the number of particles in a 1 x 1 x A cell, 
(1 — h)X = A — a, will have to be an integer. However, be- 
cause a ~ 2/3 in the limit of small doping, the condition 
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FIG. 6: (Color online) Charge density p (top) and order pa- 
rameter m (bottom) vs. U. The system has doping h = 0.05. 
Each curve is a ID cut along z-direction, the direction of the 
modulation. Beyond U c , the (100)-SDW/CDW amplitudes 
increase with U and the solution evolves from a sinusoidal 
wave to domain walls. The CDW amplitude is much weaker 
than that of the SDW. 



cannot be satisfied, and the system is necessarily metal- 
lic. A related way to see this is to consider the case of 
domain wall states, for example U = 3.5 in Fig.[6| Inside 
each domain wall (nodal region) are localized holes whose 
integrated (along the direction of the modulation) den- 
sity is a. Thus the domain wall as a whole will act as a 
quasi-2D liquid of holes with non-integer density. We will 
discuss the corresponding momentum space signature in 
the next section. 



B. SDW modulation along the (lll)-direction 



As shown above and further discussed in Sec. |IVC[ an 
orientation of the SDW other than (100) is not the so- 
lution in the proximity of U c . However, when the inter- 
action grows larger, other Fermi liquid instabilities be- 
come possible. This fact is clearly displayed when cal- 
culations on supercells commensurate with the optimal 
(100) -wavelength for a given U do not yield a state with 
(100)-SDW order. Figure [7] shows the occurrence of such 
a case in a calculation with ft = 1/8 and U = 5.0, for 
which A(ioo) = 5.5. The 16 x 16 x 22 supercell should 
have precisely accommodated 4 nodal planes of the or- 
der parameter, but rather than doing so, the random 
search produces the lower energy solution shown in the 
left panel of the figure. 

To search for the solution at higher J7, we investigate 
unidirectional SDW's with modulation lying along either 
the (110)- or (lll)-direction. An example is given in 
Fig. [8j The energies from constrained searches using the 
second approach are shown as a function of A for a scan 
of U values. It is seen that, at and above U = 4.5, the 
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FIG. 7: (Color online) Order parameter m at h = 1/8 and 
U = 5.0 in a 16 x 16 x 22 supercell (left) and a 16 x 16 x 16 
supercell (right). Though the supercell of 16 x 16 x 22 is 
commensurate with the optimal (100) -wavelength, the ran- 
dom search produces a lower energy solution. The minimum 
energy solution is an SDW along the (11 Indirection, which 
is correctly reproduced by a random search in a 16 x 16 x 16 
supercell as shown on the right. 
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FIG. 8: (Color online) Energies per site of SDW in (100)-, 
(110)- and (lll)-directions vs. A for h = 1/8. At and above 
U — 4.5, the lowest energy state is modulated along the (111)- 
instead of (lOO)-direction. 



lowest energy state is given by a (lll)-SDW, instead of 
the (100)-order at lower U. For U = 5, the minimum en- 
ergy solution is correctly reproduced by a random search 
in a 16 x 16 x 16 supercell as shown in the right panel of 
Fig. [7[ In our searches, (llO)-direction SDW's are never 
found to be the global ground state. 

By repeating the same procedure, we construct the 
equation of states for U = 3 and U = 4 contained in 
Fig. [9j At U = 3, there is no density regime where the 
(lll)-SDW is the global ground state. In contrast, for 
U = 4, a discontinuous transition from (100) to (111) 
occurs around n = 0.9 with a small coexistence region. 
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FIG. 9: (Color online) Ground-state energy per site from 
constrained search of (100)-, (110)- and (lll)-SDW at U = 3 
(left) and U — 4 (right). A linear common shift has been 
applied to the energies to highlight the convexity and the 
different trends. At U = 3, (lll)-SDW is not the global 
ground state. (The inset shows a zoom of the energy difference 
between (111)- and (100)-SDW states as n ^ 1.) At U = 4, 
the ground state is (lll)-SDW for n > 0.92. 



In both cases the low-doping ground state is character- 
ized by a linear energy-density dispersion. This bears 
two important consequences. First, contrary to what is 
observed using variational states with a uniform spiral or- 
der parameter 25 , there is no sign of phase separation into 
a half-filled, AFM region and a hole-rich region. Second, 
the effective interaction between domain walls is short- 
ranged and their precise location in the hole-diluted limit 
is therefore irrelevant as long as they stay sufficiently far 
apart. 

We find OL{m) = 1 at any density for which the (111)- 
SDW is the ground state. The (lll)-SDW states are 
fully gapped, owing to the integer value of a/m), in con- 
trast to the metallic behavior of the (100) states. Upon 
increase of U at a constant h, the structural transition 
is therefore always accompanied by a metal-to-insulator 
transition. We have verified, for selected cases, that ran- 
dom searches on larger supercells with sizes commensu- 
rate to the optimal wave- vector always return unidirec- 
tional SDW's with the same predicted wavelength and 
orientation. This provides a strong indication that the 
character of intermediate U instabilities remains that of 
a unidirectional SDW. Thus, as we increase U at con- 
stant density, the system is always expected to undergo 
a discontinuous transition from a (100)- to a (lll)-SDW 
ground state. 



C. A variational pairing ansatz 

The pairing model is a variational ansatz that has 
proved extremely helpful in rationalizing the properties 
of SDW 's in the mean-field treatment of the electron 
ga J 22 | 2e | and the 2D Hubbard modePP. Similarly here, 



the model helps to explain the numerical results and pro- 
vides a simple conceptual framework that captures the 
essential physics of the SDW states in 3D and in the 
crossover regime discussed in the next section. We first 
summarize the formalism, and then apply it to the case 
of (100)-SDW at modest U, followed by the (lll)-SDW. 

At low U and small /i, the pairing model is defined by 
spin-orbitals of the form 



(12) 



The construction requires excitations to outside the 
Fermi sea. It is reminiscent of the ansatz used to con- 
struct the BCS pairing states for attractive interactions, 
except that here the tendency for small excitations is 
perhaps more "natural", because of the repulsive inter- 
action. A collection of spin pairs in orbitals given by 
Eq. (12) leads to a uniform charge density, p(R) = n, 
and a spin density of the form 



8(B) 



N ^ 

ken 



ak cos(qk • R) 



(13) 



with ak = ^k^k- The region 1Z over which k is summed 
will be closely related to the non-interacting Fermi sea, 
and preserving the volume of 47r 3 n, but will in general be 
slightly modified from the variational optimization, as we 
further discuss below. To ensure orthogonality amongst 
the spin-orbitals, qk must be such that k + qk 7Z and 
k + 7^ k' + qk' . The corresponding potential energy 
per site is then given by 



V 



u_ 

47V ■ 



Un 2 -^^ 5 2 (R). 

R 



(14) 



The potential energy lowering relative to the paramag- 
netic (PM) solution is thus: 



AV 



U 
"iV 2 



^ a k a k ' [<S(qjc + q k ') + <S(q k - q k ')] > 

k,k'G^ 

(15) 



where the Kronecker S is intended as periodic on the 
reciprocal lattice, i.e., modulo 2tt in any direction. 



Equation ( 15 ) makes it clear that the maximum reduc- 
tion in V is achieved by having as many pairs as possible 
with qk's which are parallel or anti-parallel to each other. 
Noting that the vector Q is perfectly nested when h = 0, 
let us consider the following explicit construction for 1Z: 
displace the half-filled Fermi surface in each octant of the 
first BZ by ±Aq/2, choosing the direction that shrinks 
the Fermi sea and with a length Aq such that the en- 
closed volume is reduced to 47r 3 n. This construction is 
illustrated in the top panel in Fig. [To] for t he case of 
the (100)-SDW (discussed next in Sec.|IVCl|). The sur- 



face of 1Z can now anchor spin pairs with one common 
pairing vector, by making ^k less than 1 in a small layer 
immediately inside the surface of 7Z, and correspondingly 

^k = V 1 -Vk| 2 > o. 
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estimate of a (Eq. ( [TT] )) 



FIG. 10: (Color online) Illustration of the pairing model for 
the (100)-SDW state in 3D. The schematic diagram is drawn 
on the actual momentum distribution from the exact numer- 
ical solution for h = 1/8 and U = 4.0. The top panel shows 
the pairing construction on the contour plot of the (110) cut. 
The white dashed lines represent the half- filled surface, across 
which is the nesting vector Q. The reconstructed surface, 
shown as magenta solid lines, is obtained by displacing the 
half- filled one along the ^-direction by a distance Ag(ioo)/2, 
as given by Eq. ( |16[ ). The nesting vector across the shifted 
surface, q(ioo) , is shown by the long solid line with arrow. The 
bottom panel shows n(k) of the k z = tt plane. This is in a 
region where the Fermi surface survives the onset of order and 
where it differs more severely from the pairing construction. 
The actual Fermi surface, seen distinctly inside the recon- 
structed surface, differs little from the non-interacting Fermi 
surface (black solid line), rik drops sharply, and no pairing is 
present in this region. 



For small /i, we can determine Aq directly from the 
construction: 



■S 



n 



/ eAq • dS = h- 



BZ 



2 ' 



(16) 



where S is the half- filled surface, eA q is the direction 
of Aq, and ^bz = (27r) 3 is the volume of the first BZ. 
Equation ( 16) implies a linear relationship between h and 
Aq and, using A = 7r/(Aq • e), provides the following 



47r 2 eAq • e 



eAq 



dS. 



(17) 



where e is a relevant Cartesian unit vector. 

Different directions of Aq lead to different reconstruc- 
tions of the non-interacting doped Fermi surface. The 
kinetic energy cost of the pairing ansatz can therefore 
be thought of as the combination of two contributions: 
the reconstruction energy due to using 1Z rather than 
the true non-interacting Fermi sea, and the kinetic en- 
ergy change due to moving particles from k (inside 7Z) to 
k + qk (outside), ie., the non-zero v^s. It is easy to see 
that, similar to 2EPHl, at sufficiently large U the potential 
energy lowering will overtake the kinetic energy increase 
for the constructions discussed here. The correct state 
is determined by maximizing the gain in the potential 
energy from pairing (larger areas near the Fermi surface 
participating with parallel Aq) while minimizing the ki- 
netic energy cost. 

The ansatz gives a clear picture for the onset of the in- 
stability. First, by its form, the model captures how the 
energetically costly CDW can be suppressed compared 
to the SDW. Second, it indicates that amongst different 
possible choices of qk, the one involving only parallel and 
anti-parallel vectors are optimal. Third, the direction of 
±Aq must be such as to lead to the minimal possible re- 
construction of the doped Fermi surface. These findings 
are in qualitative agreement with the numerical results 
obtained by the solution of the SCF equations: 1) the 
SDW is much stronger than the accompanying CDW, 
2) the Fermi surface reconstructs in a way to enhance 
pairing and the SDW order tends to be unidirectional as 
a result, 3) more drastic reconstructions are only possi- 
ble with larger U . Much quantitative information can 
be obtained with straightforward calculations using this 
model, as we discuss next for the (100)- and (lll)-SDW 
states, respectively. 



1. Analysis of the (100) -SDW 

Among all directions, only a Aq along the (100)- 
direction causes the Fermi surface in each octant to 
shrink equally and this, it can be shown, leads to the 
minimum kinetic energy increase at small h. An SDW 
with (100) modulation is thus the lowest energy solution 
at low h and just above U c , consistent with the results 
from explicit solutions of the SCF equations in Sec. IV A| 
We numerically calculate the projected area along the 
(lOO)-direction (shown in the right middle panel in Fig. [I]) 



and obtain a 



(100) 



0.63, in very good agreement with 



the exact results from direct solutions shown in Fig. [5j 
where a^ 10 o) — 0.66. That the estimated value is slightly 
smaller is consistent with the presence of surviving Fermi 
surface inside the reconstructed doped Fermi surface as 
seen in Fig. [lOj 
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FIG. 11: (Color online) Contour plots of the single-particle 
spectral function evaluated at the Fermi energy; on the (110) 
cut (top) and on the k z = 0,7r/2,7r planes (bottom from left 
to right) for a system with U — 2.7 and h = 1/8. A large part 
of the Fermi surface survives, except for areas around the hot 
spots where it is most energetically favorable for pairing. 



A direct comparison between the pairing model and 
the exact SCF solution can also be made in momentum 
space. We will identify the Fermi surface in the numeri- 
cal solution from mean-field theory as the locus of points 
where = 0.5. Depending on the system and the value 
of U , the momentum distribution of the exact mean- field 
ground state can maintain a true Fermi surface, char- 
acterized by a discontinuity in nk, or have it smeared 
out by large pairing amplitudes (i.e., close to l/y/2 
near the boundary of 1Z) . The two scenarios can occur in 
the same system at different k values. The identification 
using = 0.5 is consistent with both. 



Figure [10| shows, in particular, that the portion of the 
Fermi surface where pairing takes place is in very good 
agreement with the construction based on the pairing 
model, which indicates that the ansatz captures the dom- 
inant ingredient of the physics of the SDW state. The 
figure also provides a direct explanation for the survival 
of the Fermi surface around k z = tt as it is there that 
the pairing construction shows large discrepancy with the 
true Fermi surface. This, in turn, implies that pairing in 
that region would be associated with too large a kinetic 
energy cost to be favorable. The absence of any gap from 
pairing at the Fermi surface in the k z = tt plane is con- 
sistent with the picture discussed earlier of a quasi-2D 
liquid within each domain wall. These effects are ampli- 
fied at smaller U : s as shown in Fig. [TT] where a larger 



■ center-point 

■ entire surface 




0.05 0.1 0.15 0.2 

h 



FIG. 12: (Color online) (a) Half-filled Fermi surface for the 
3D Hubbard model in one of the octants (left) and (b) aver- 
aged Ag( 10 o) over different areas vs. doping (right). The blue 
solid line in (b) is calculated at (— 7r/2, — 7r/2,7r/2), shown as 
the blue dot in (a). The dashed/dotted lines in (b) are aver- 
aged over varying areas as indicated by the circles with the 
same color/style in (a). The black line in (b) shows the aver- 
aged value over the entire surface, which leads to the estimate 



of a 



(ioo) 



0.63 discussed in the text. 



part of the Fermi surface survives the onset of order; the 
parts that do not survive are in areas around the "hot 
spots" k ~ (±7r/2, ±7r/2, ±7r/2), where the pairing con- 
struction and the doped Fermi surface are most similar, 
and the change in is at a minimum implying a closer 
proximity to a perfect common paring vector. 

These understandings allow quantitative explanations 
of all the features of the data on a (inn) shown in Fig. [5j 
For example, Fig. [12] shows that the distance along z- 
direction between the doped and half-filled Fermi surfaces 
is at a minimum around the hot spots. Given that such 
distance equals Ag/2 in the pairing construction, one 
finds that the local Aq value is smaller at the hot spots 
than when computed as the average distance over the en- 
tire surface. Hence, when only the hot spots are involved 
in pairing, a larger a (mo) results, as seen at lower values 
of U just above U c . Obviously, the pairing reconstruction 
becomes increasingly accurate as h approaches 0. In this 
limit, one therefore finds the increasingly smaller U c and 
the faster convergence (in U) to the saturated value of 
a <ioo) ~ 2/3 shown in Fig. [5] 



2. Application to the {111) -SDW 

By taking a displacement Aq along the (11 Indirection, 
we can apply the pairing construction straightforwardly 
to the diagonal-modulated SDW. Analogously to the 
(100) case, Fig. 13 shows a remarkable agreement be- 



tween the shifted half-filled surface and the calculated 
Fermi surface of the SDW. The two cuts in the figure 
clearly show broken cubic symmetry, where the Fermi 
surface in one pair of the octants is further away from 
the half-filled one so as to share the common modula- 
tion wave-vector Aq/m) with the other three pairs. The 
(lll)-SDW state offers, in this respect, a particularly 
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FIG. 13: (Color online) Schematic illustration of the pairing 
model for the (lll)-SDW state, shown on the contour plots of 
the (110) (top) and the (110) (bottom) cuts of nkt for a sys- 
tem of h = 1/8 and U = 5.0. The white dashed lines are the 
half- filled surface, across which is the nesting vector Q. The 
magenta solid lines show the pairing construction, obtained 
by shifting the half-filled surface along [Ill]-direction by a 
distance of A^m)/2. In the upper panel, q(in) and q/m) 
give two equivalent representations (differing by a reciprocal 
lattice vector) of the pairing vector across the reconstructed 
Fermi surface. Note the asymmetry between the two diagonal 
directions in the upper panel. 



clear example where Fermi surface reconstruction can be 
observed. It also shows how an accurate experimental 
characterization of the momentum distribution in opti- 
cal lattices can be used to characterize the band structure 
and pairing at the Fermi surface which, in turn, provides 
momentum space evidence on the real space character of 
the SDW. 



Using Eq. (17) we have estimated the wavelength of 
the (lll)-SDW, an d foun d a {111) = 0.93. The exact SCF 
calculations in Sec. IV B showed, instead, that a/m) is 
precisely pinned at 1 in a fairly large regime of U. This 
somewhat large discrepancy is a consequence of the nat- 
ural tendency of the system to "lock" the integrated den- 
sity of holes per nodal region to 1 whenever the topology 
of the non-interacting doped Fermi surface is such that 



this is not energetically too costly. When doing so the 
system benefits from both pairing and band-insulating 
mechanisms to gap the entire Fermi surface and further 
lower the ground-state energy. 



V. DIMENSIONAL CROSSOVER RESULTS 

A. Results from full numerical HF solutions 

The mean-field ground state of the doped 2D Hubbard 
model shares many similarities with its 3D counterpart. 
Just above U c , the 2D system develops a sinusoidal SDW 
with a modulating wave along the (lO)-direction and a 
much weaker accompanying CDW. As U is increased, 
the SDW increases its amplitude before the SDW state 
eventually changes into a collection of weakly interacting 
domain walls. Above a certain U, there is a discontinuous 
transition to a phase where the modulation is along the 
(1 Indirection. The crossover from SDW to domain walls 
occurs before the (10) to (11) transition at small h, but 
after at larger /PH. A peculiarity of the 2D case, due to 
the special topology of the 2D half-filled surface, is that 
a = 1 and the system is an insulator regardless of doping, 
U or direction of the modulation wave- vector apart from 
a region close to U c . 

By controlling the distance between square lattice lay- 
ers, optical lattice experiments allow the study of the evo- 
lution of the system as it crosses over from 2D to 3D. This 
situation is theoretically described by an increase of t± in 
Hamiltonian ([!]) and the question, within mean-field the- 
ory, concerns the ensuing evolution of the ground-state 
properties. The pairing model and the arguments de- 
scribed in Sec. IV. C remain valid in the crossover regime. 
We thus restrict our investigation to unidirectional SDW 
ground states, although we did use the first approach 
to carry out some searches, finding no additional struc- 
tures. As in 3D, we verify that the SDW solution with 
minimum energy, identified using the second approach, 
can be obtained by the first approach in a large super- 
cell that is commensurate with the optimal wavelengths, 
even when starting from random initial guesses. SDW in 
directions different from (100) or (111) are not found to 
be the global ground state for any value of t±. 

Results are summarized in the t±-U mean-field phase 
An overall increase in the critical U 



diagram of Fig. 14 



values is seen as doping increases, as a result of a greater 
deformation from the perfectly nested half-filled Fermi 
surface and the need for more excitations to achieve re- 
construction of the Fermi surface for pairing. As before, 
numerical calculations focus on small doping (h ^ 0.2) 
and low to intermediate interactions (U ^ 5.5), where 
mean-field theory can be expected to be more accurate. 
Upon increase of [/, and similarly to 3D, the system un- 
dergoes a first transition to a (100)-SDW state followed 
by a second, discontinuous transition to a (lll)-SDW 
state. The absence of cubic symmetry away from t± = 1 
causes the modulation wave-vector for (100)-SDW to lie 
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FIG. 14: (Color online) Mean-field phase diagram of the 
ground state in the crossover regime. Phase boundaries 
for several values of doping are indicated by symbols. The 
lines are to guide the eye. Solid lines separate the PM 
phase from the AFM phase, and dashed lines show transi- 
tions from (100)- to (lll)-SDW. The inset plots the value of 

^(100)^(111) /^PM^(IOO) • 



in the x?/-plane. This is because the elongation of the 
Fermi surface along z-direction (as illustrated in Fig. [I]) , 
meaning that Aq along z-direction leads to less surface 
area for pairing than along x- or ^-directions. Wave- 
vectors along the (11 Indirection continue, on the other 
hand, to remain equivalent under the symmetry opera- 
tion of the tetragonal group. 

The critical value of the interaction strength for the 
transition from the paramagnetic (PM) phase to (100)- 
SDW, £7pm^(ioo)? monotonically increases from 2D to 
3D, due to the wider band width and smaller density of 
states at the Fermi energy for larger t±. The transition 
values decrease to when h approaches as U c = for 
the half- filled system at any t±. The critical U value for 
the transition from (100)- to (lll)-SDW, E/<ioo)-Kiii)> 
has a lower bound lying close to the h = 1/32 line in the 
figure, so that no (lll)-SDW exists below U — 3 regard- 
less of the smallness of h and the value of t± . In contrast 



with Upm- 



-(100)' 



100)^(111) 



displays a non-monotonic 



behavior with t± whose origin we will address in the next 
section. 

The evolution of the modulation wavelength is sum- 



marized in Fig. 15 in terms of a as a function of t±. 
Numerical results are obtained at U values around the 
transition line in Fig. 



14 



^ e.g., just below U(ioo>-+<iii) for 
(100). In a small winclow around t± = 0, a(ioo) remains 
1, as a direct consequence of the insulating character of 
the 2D solution. Apart from this region, there is a grad- 
ual decrease of a(ioo) with t± (discussed further in the 
next section) and a weakening of the intensity of the or- 
der parameter (bottom panel of Fig. 15). As we have 
already remarked in the 3D results, the pinning of a at 1 
in the (111) phase, and the corresponding insulating be- 
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FIG. 15: (Color online) Top panel: dependence of the modu- 
lation wavelength on t± . Numerical results for a (data points 
with error bars) at various dopings are compared with theo- 
retical estimates (solid lines). Bottom panel: evolution of the 
order parameter on a 48 x 10 x 10 supercell for h = 3/25, 
U — 3.5 as a function of t±. As t± increases, a<ioo) remains 
at 1 in a very small region close to 2D, before decreasing 
monotonically, while CK(m) is locked at 1. 



havior, are a result of a further stabilization of the ground 
state caused by coexisting magnetic and band-insulating 
effects. 



B. Pairing model discussion 



The dependence of a on t± is captured by Eq. (17), 
as the projection of the Fermi surface S(tj_) along Aq. 
The middle and bottom rows of Fig. [l] show how the 
projected surfaces, along the (100)- and (11 Indirections 
respectively, shrink as t± is increased. Results from ex- 
plicit calculations using Eq. ( fl7| ) are shown as continuous 
lines in Fig. [15] They display a correct trend but a con- 
sistently underestimation of wavelengths. The quantita- 
tive disagreement is not surprising. The most significant 
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reason behind it is the tendency of a to be locked at 1, 
on which we have already commented and which appar- 
ently involves more global considerations than contained 
in the pairing model. The smaller discrepancy for the 
(lOO)-SDW outside the immediate vicinity of 2D, which 
increases for higher h, is due to the surviving FS that 
remains inside the reconstructed doped FS, as we have 
already remarked in the 3D results. 

We next address the origin of the non-monotonic be- 
havior of £7(100)^(111)- This is the result of two compet- 
ing factors. On the one hand, the increasing band width 
with dimensionality leads to an increase of U c , as demon- 



strated in the monotonicity of C/pm- 



-(100) 



. On the other 



hand, the geometrical properties of the Fermi surface are 
such that the angle between Aq/m\ and the Fermi sur- 
face in some of the octants is small when t± is small. 
This means that the displaced Fermi surfaces in those 
octants will remain close to the half-filling counterparts 
in the construction of 1Z (Sec. IV C). As a result, the 
other components of the Fermi surface must be displaced 
further to preserve total volume, causing a more uneven 
reconstruction which requires more excitations from the 
non-interacting Fermi sea, hence larger U c . The 2D case 
offers an extreme example of this as it is characterized 
by a large fraction of the reconstructed doped surface re- 
maining exactly pinned on the half-filled one 20 . To sepa- 
rate this factor from that of the band width, we examine 



(100) 



Fig. 14 



;iii)/E/pm->>(ioo)> which is plotted in the inset of 
A monotonic decrease is seen with t±. There- 
fore, £7(ioo) ->> (in) fi rs t decreases and then increases as t± 
goes from to 1. 



VI. SUMMARY AND DISCUSSION 

This work addressed quantitative aspects of possi- 
ble inhomogeneous magnetic phases of the 3D Hubbard 
model that emerge as the average density deviates from 
one particle per site. Because of the ease with which ex- 
periments are expected to be able to transition between 
the 2D and 3D regimes, we also studied the evolution 
of the inhomogeneous ground state as a function of the 
hybridization between parallel layers of square lattices. 
Within mean-field theory, we have shown that the lead- 
ing instability of the PM ground state is an SDW with 
long wavelength modulation along the (100) -direction. 
No tendency toward phase separation was seen, even at 
small values of doping. The system remains metallic, re- 
gardless of the proximity to half-filling, because of a non- 
integer density of holes per wavelength of modulation. 
This density is largely determined by an entirely geo- 
metric property of the Fermi surface: its projected area 
along the direction of modulation. At larger U values, the 



ground state continues to be a unidirectional SDW, but 
with (lll)-orientation. This phase is insulating and char- 
acterized by a significant distortion of the momentum 
distribution. Such distortion leads, quite naturally, to 
the identification of a reconstructed Fermi surface whose 
observation in optical lattice experiments should be fea- 
sible. 

We showed that much of these results can be under- 
stood by a simple variational ansatz with pairing orbitals 
formed by a linear combination of two plane waves. By 
placing a pair of up- and down-spin particles into a pair 
of such orbitals, an SDW is formed with constant charge 
density. Straightforward analysis of the energetics from 
this ansatz leads to quantitative predictions of the wave- 
length and nature of the SDW modulation which are ver- 
ified by our direct numerical solutions of the SCF equa- 
tions. 

The true many-body ground state will modify the 
mean-field solutions in several ways. For example, quan- 
tum Monte Carlo calculations in periodic simulation cells 
will restore translational invariance, and the inhomo- 
geneities seen here will be manifested in spin-spin correla- 
tions and such. A deeper issue is the possible existence of 
additional competing instabilities once a fuller treatment 
of quantum fluctuations is included. Certainly, the ten- 
dency for magnetic inhomogeneous order is exaggerated 
in mean-field theory, and a more accurate description of 
the many-body correlation at a certain U value tends 
to be given by the mean-field results at a significantly 
weaker U. However, as we have shown in 2D, mean- 
field theory appears to capture the correct basic picture 
of the magnetic correlations when compared to quantum 
Monte Carlo results^^oj. This indicates that the results 
in the present paper can provide a useful framework for 
understanding the magnetic correlations in 3D and in 
the crossover regime for weak to intermediate interaction 
strengths. 

Apart from the obvious omissions inherent in the 
mean-field approximation, this study has not addressed 
the fact that experiments are performed in the presence 
of a confining potential. Nor have we addressed how the 
situation is modified by a finite magnetization. Gener- 
alization of the present approach to address such issues 
will be valuable, and technically straightforward. 
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